function [dxdt,s]=rhs(t,states,p)
    s=struct();
    x1=states(1); y=x1;
    x2=states(2);
    dxdt=zeros(size(states));

    if t<p.Ts
        mu=t/p.eps/(1-(p.rho)*t+p.rho*p.Ts);
    else
        mu=1/p.eps;
    end
    yd=p.yd(t);
    z1=mu*(y-yd);
    rho1hat=states(3);
    chi=1/(1-z1^2);
    theta1=zeros(2,1);
    phi1=zeros(2,1); 
    xi1=[phi1;1];
    sigma1=p.sigma1(t);
    sigma2=p.sigma2(t);
    eta11=chi*(y-yd)^3/sqrt((z1*chi*(y-yd))^2+sigma1^2);
    eta12=z1*chi*mu*(xi1'*xi1)/sqrt(z1^2*chi^2*mu^2*(xi1'*xi1)+sigma1^2);
    eta1=eta11+eta12;
    alphabar1=p.k1*(y-yd)*chi+rho1hat*eta1;
    alpha1=-z1*chi*mu*alphabar1^2/sqrt((z1*chi*mu*alphabar1)^2+sigma1^2); %(18)
    dxdt(3)=p.gamma1*z1*chi*mu*eta1-p.gamma1*p.sigmarho1(t)*rho1hat;   %(19)
    
    
    % virtual adaptive control laws (18)-(19) and (34)-(35)


    % the adaptive controller (41)-(42)
    z2=x2-alpha1;
    ds=diffs([t y yd mu sigma1 rho1hat]);
    l2=-ds.alpha1_rhohat1*dxdt(3);

    M=p.M(t);K=p.K(t);C=p.C(t);
    kappa=p.kappa(t);
    ubar=p.ubar(t);
    g1=0;
    g2=1/M;
    theta2=[-K/M;-C/M];
    phi2=[x1;x2];

    Theta21=[theta2;theta1;g1];
    xi21=[phi2;0;-ds.alpha1_x1*x2];
    xi22=[1;-ds.alpha1_x1];
    xi23=[-ds.alpha1_yd;-ds.alpha1_mu;-ds.alpha1_sigma1;];
    xi2=[xi21;xi22;xi23];
    xibar2=[xi2;1;z1;l2];
    eta2=z2*(xibar2'*xibar2)/sqrt(z2^2*(xibar2'*xibar2)+sigma2^2);
    rhohat2=states(4);dxdt(4)=p.gamma2*z2*eta2-p.gamma2*p.sigmarho2(t)*rhohat2;
    alphabar2=p.k2*z2+rhohat2*eta2;
    nu=-z2*alphabar2^2/sqrt((z2*alphabar2)^2+sigma2^2);

    % system 51
    % the planet (1)
    % unknown actuator failure (2)
    
    dxdt(1)=x2+p.d1(t);
    dxdt(2)=g2*(kappa*nu+ubar)+theta2'*phi2+p.d2(t);
    s.yd=yd;
end